function [varargout] = alphasim(ddim, opts)
% alphasim  - simulate noise data to estimate cluster threshold
%
% FORMAT:       [at = ] alphasim(ddim [, opts])
%
% Input fields:
%
%       ddim        data dimension (1x3 integer values)
%       opts        optional settings
%       .clconn     connectivity of clusters, ('face', {'edge'}, 'vertex')
%       .conj       conjunction simulation (1x1 double, number of maps)
%       .fftconv    boolean flag, use FFT convolution (default: false)
%       .fwhm       FWHM kernel sizes (default: [2, 2, 2])
%       .mask       boolean mask (size must be == ddim!, default: none)
%       .niter      number of iterations, default: 1000
%       .pbar       either BVQXprogress or BVQXfigure:XProgress object
%       .thr        applied (raw) threshold(s), default: p<0.001
%
% Output fields:
%
%       at          optional output table
%
% Note: other than AFNI's AlphaSim, the data is considered to be
%       iso-voxel for the default kernel, but that can be altered
%       accordingly by changing the kernel!

% Version:  v0.8a
% Build:    9102220
% Date:     Oct-22 2009, 8:43 PM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
   ~isa(ddim, 'double') || ...
    numel(ddim) ~= 3 || ...
    any(isinf(ddim) | isnan(ddim) | ddim < 1 | ddim > 256)
    error( ...
        'BVQXtools:BadArgument', ...
        'Missing or invalid ddim argument.' ...
    );
else
    ddim = round(ddim);
end
if nargin < 2 || ...
    isempty(opts)
    opts = struct;
elseif ~isstruct(opts) || ...
    numel(opts) ~= 1
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid opts argument.' ...
    );
end
if ~isfield(opts, 'clconn')
    opts.clconn = 'edge';
elseif ~ischar(opts.clconn) || ...
   ~any(strcmpi(opts.clconn(:)', {'edge', 'face', 'vertex'}))
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid opts.clconn field.' ...
    );
else
    opts.clconn = lower(opts.clconn(:)');
end
switch (opts.clconn(1))
    case {'e'}
        clconn = 2;
    case {'f'}
        clconn = 1;
    otherwise
        clconn = 3;
end
if ~isfield(opts, 'conj')
    nconj = 1;
elseif ~isa(opts.conj, 'double') || ...
    numel(opts.conj) ~= 1 || ...
    isinf(opts.conj) || ...
    isnan(opts.conj) || ...
    opts.conj < 1 || ...
    opts.conj > 5
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid opts.conj field.' ...
    );
else
    nconj = floor(opts.conj);
end
if ~isfield(opts, 'fftconv')
    opts.fftconv = false;
elseif ~islogical(opts.fftconv) || ...
    numel(opts.fftconv) ~= 1
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid opts.fftconv field.' ...
    );
end
fftconv = opts.fftconv;
if ~isfield(opts, 'fwhm')
    opts.fwhm = [2, 2, 2];
elseif ~isa(opts.fwhm, 'double') || ...
    numel(opts.fwhm) ~= 3 || ...
    any(isinf(opts.fwhm) | isnan(opts.fwhm) | opts.fwhm <= 0 | opts.fwhm(:)' > ddim)
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid opts.fwhm field.' ...
    );
else
    opts.fwhm = opts.fwhm(:)';
end
kcell = { ...
    smoothkern(opts.fwhm(1), 0), ...
    smoothkern(opts.fwhm(2), 0), ...
    smoothkern(opts.fwhm(3), 0)};
kern = {zeros(numel(kcell{1}), numel(kcell{2}), numel(kcell{3}))};
kern = kern(1, [1, 1, 1]);
kern{1}(:, (numel(kcell{2}) + 1) / 2, (numel(kcell{3}) + 1) / 2) = kcell{1};
kern{2}((numel(kcell{1}) + 1) / 2, :, (numel(kcell{3}) + 1) / 2) = kcell{2};
kern{3}((numel(kcell{2}) + 1) / 2, (numel(kcell{2}) + 1) / 2, :) = kcell{3};
if ~isfield(opts, 'mask')
    opts.mask = ([] > 0);
elseif ~islogical(opts.mask) || ...
   ~isequal(size(opts.mask), ddim)
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid opts.mask field.' ...
    );
end
mask = opts.mask;
if ~isfield(opts, 'niter')
    niter = 1000;
elseif ~isa(opts.niter, 'double') || ...
    numel(opts.niter) ~= 1 || ...
    isinf(opts.niter) || ...
    isnan(opts.niter) || ...
    opts.niter < 1 || ...
    opts.niter > 1e6
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid opts.niter field.' ...
    );
else
    niter = round(opts.niter);
end
if ~isfield(opts, 'pbar') || ...
    numel(opts.pbar) ~= 1 || ...
   ~any(strcmpi(class(opts.pbar), {'BVQXfigure', 'BVQXprogress'}))
    opts.pbar = [];
end
if ~isfield(opts, 'thr')
    thr = 0.001;
elseif ~isa(opts.thr, 'double') || ...
    isempty(opts.thr) || ...
    numel(opts.thr) > 16 || ...
    any(isinf(opts.thr(:)) | isnan(opts.thr(:)) | opts.thr(:) <= 0 | opts.thr(:) > 0.5)
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid opts.thr field.' ...
    );
else
    thr = opts.thr(:)';
end
nthr = numel(thr);

% compute z-threshold
if isempty(mask)
    pdim = prod(ddim);
else
    pdim = sum(opts.mask(:));
end
rpos = ceil(pdim * thr);
zthr = -custom_tinv(thr, 1e7);

% create counting arrays
cc = zeros(nthr, 1000);
fc = zeros(nthr, niter);

% prepare convolution FFT kernel if required
if fftconv
    kern = smoothkern(opts.fwhm, 0);
    kdim = size(kern);
    kdimh = floor(kdim ./ 2);
    fftkern = zeros(ddim);
    ddh = round((ddim + 1) / 2);
    fftkern(ddh(1)-kdimh(1):ddh(1)+kdimh(1), ...
            ddh(2)-kdimh(2):ddh(2)+kdimh(2), ...
            ddh(3)-kdimh(3):ddh(3)+kdimh(3)) = kern;
    fftkern = fftn(fftkern);
end

% test BVQXprogress
if niter >= 50
    if isempty(opts.pbar)
        try
            pbar = BVQXprogress;
            BVQXprogress(pbar, 'setposition', [80, 200, 640, 36]);
            BVQXprogress(pbar, 'settitle', 'Running alphasim...');
            BVQXprogress(pbar, 0, sprintf('0/%d iterations, %d thresholds...', ...
                niter, nthr), 'visible', 0, 1);
            pbarn = '';
            pst = niter / 100;
            psn = pst;
        catch
            pbar = [];
            psn = Inf;
        end
    else
        pst = ceil(niter / 500);
        psn = pst;
        pbar = opts.pbar;
        pbar.Progress(0, sprintf('alphasim: 0/%d iterations, %d thresholds...', niter, nthr));
        pbarn = 'alphasim: ';
    end
else
    pbar = [];
    psn = Inf;
end

% extend ddim
ddim(4) = nconj;

% run loop
for n = 1:niter
    
    % convolution type
    if ~fftconv
        if nconj == 1
            r = conv3d(conv3d(conv3d(randn(ddim), kern{1}), kern{2}), kern{3});
        else
            r = randn(ddim);
            for nc = 1:nconj
                r(:, :, :, nc) = conv3d(conv3d(conv3d(r(:, :, :, nc), ...
                    kern{1}), kern{2}), kern{3});
            end
        end
    else
        if nconj == 1
            rf = fftn(randn(ddim));
            rf = rf .* fftkern;
            r = fftshift(ifftn(rf));
        else
            r = randn(ddim);
            for nc = 1:nconj
                rf = fftn(r(:, :, :, nc));
                rf = rf .* fftkern;
                r(:, :, :, nc) = fftshift(ifftn(rf));
            end
        end
    end
    
    % do for each threshold
    for tc = 1:nthr
    
        % no conjunction
        if nconj == 1
            
            % scale so that thr
            if isempty(mask)
                rs = sort(r(:));
            else
                rs = sort(r(mask));
            end
            r = r .* (zthr(tc) / rs(rpos(tc)));

            % mask
            if ~isempty(mask)
                rf = r .* mask;
            else
                rf = r;
            end

        % conjunctions
        else
            
            % scale each map
            for nc = 1:nconj
                
                % scale so that thr
                rp = r(:, :, :, nc);
                if isempty(mask)
                    rs = sort(rp(:));
                else
                    rs = sort(rp(mask));
                end
                rp = rp .* (zthr(tc) / rs(rpos(tc)));

                % mask
                if ~isempty(mask)
                    if nc == 1
                        rf = rp .* mask;
                    else
                        rf = sign(rf) .* (sign(rf) == sign(rp)) .* ...
                            min(abs(rf), abs(rp .* mask));
                    end
                else
                    if nc == 1
                        rf = rp;
                    else
                        rf = sign(rf) .* (sign(rf) == sign(rp)) .* ...
                            min(abs(rf), abs(rp));
                    end
                end
            end
        end

        % cluster frequency
        cf = clustercoordsc(rf >= zthr(tc), clconn);

        % largest cluster
        if ~isempty(cf)
            mc = max(cf);
        else
            mc = 0;
        end

        % extend array if necessary
        if mc > size(cc, 2)
            cc(1, mc + ceil(size(cc, 2) / 12)) = 0;
        end

        % put into frequency arrays
        fc(tc, n) = mc;
        for nc = 1:numel(cf)
            cc(tc, cf(nc)) = cc(tc, cf(nc)) + 1;
        end
    end
    
    % update progress bar
    if n >= psn && ...
       ~isempty(pbar)
        pbar.Progress(n / niter, sprintf(...
            '%s%d/%d iterations, %d thresholds...', pbarn, n, niter, nthr));
        pbar.Visible = 'on';
        psn = psn + pst;
    end
end

% close progress bar
if ~isempty(pbar) && ...
    isempty(opts.pbar)
    closebar(pbar);
end

% get size and data
mf = max(fc, [], 2);
cc = cc(:, 1:max(mf));

% prepare output
tout = cell(nthr, 1);
for tc = 1:nthr
    hf = hist(fc(tc, :), 1:mf(tc));
    hfs = cumsum(hf(:));
    hx = cc(tc, 1:mf(tc)) .* (1:mf(tc));
    hxs = [0; hx(:)];
    hxs(end) = [];
    ht = sum(hx);
    sc = sum(cc(tc, 1:mf(tc)));
    ccx = cc(tc, :)';
    tout{tc} = [(1:mf(tc))', ccx(1:mf(tc)), cumsum(ccx(1:mf(tc))) ./ sc, ...
        thr(tc) .* (sum(hx) - cumsum(hxs(:))) ./ ht, hf(1:mf(tc))', ...
        1 - ([0; hfs(1:mf(tc)-1)]) ./ niter];
end

% output variable or table
if nargout < 1
    for tc = 1:nthr
        if nthr > 1
            disp(' ');
            disp(sprintf('Uncorrected threshold: p<%f', thr(tc)));
            disp('-----------------------------------------------------------');
        end
        disp('Cl Size  Frequency  CumProbCl  p / Voxel  MaxFreq   Alpha ');
        disp(sprintf('%7d  %9d  %9.7f  %9.7f  %8d  %7.5f\n', lsqueeze(tout{tc}')));
    end
else
    if nthr == 1
        varargout{1} = tout{1};
    else
        varargout{1} = tout;
    end
end
